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I . INTRODUCTION 



There are many theories to aid the airfoil design 
process by computational methods, because of the desire to 
reduce the number and cost of wind tunnel tests. A still 
largely unresolved question is the problem of flow sepa- 
ration. Because the classical boundary layer approximation 
cannot be applied to separated flow calculations, engineers 
have tried to overcome this limitation by developing 
viscous/ inviscid interaction approaches or to develop direct 
solutions of the Navier-Stok.es equation. 

The purpose of this thesis is to demonstrate the 
capability of the viscous/ inviscid interaction method by 
applying Cebeci's interactive computer program to a single 
airfoil ( FX 63-137) at three low Reynolds numbers and by 
comparing the results with experimental data. 

Chapter II explains the boundary layer theory. The 
boundary layer equations are derived and the turbulence 
models are introduced. Also, this chapter includes the 
prediction of transition boundary layer calculations and 
flow separation. 

Chapter III introduces the interaction methods. 

Three weak, interaction methods are explained briefly and the 
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simultaneous method is presented as a strong interaction 
me thod . 

Chapter IV describes Cebeci's interactive computer 
program. Input/Output data description and JCL files are 
included. Also, the results of the application of Cebeci's 
program are discussed. 
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BOUNDARY LAYER THEORY 



II . 



A. DERIVATION OF BOUNDARY LAYER EQUATIONS 

Generally, the thickness of the boundary layer increases 
with viscosity, or it is possible to state that it decreases 
with viscosity, or it is possible to state that it decreases 
as the Reynolds number increases. From exact solutions of 
the Nav ier -Stokes equations, it was seen that the 
boundary-layer thickness (<T ) is proportional to the square 
root of kinematic viscosity ())); 

cT(x) wo / l) X 

V u 

where x is the distance from the leading edge of a flat 
plate. Using the local Reynolds number 



constant - property flow without body forces and leave the 
stresses unspecified so that the results apply to laminar or 
turbulent flow. If simplifications are to be introduced 
into Navier -Stokes equations, it must be assumed that the 
boundary layer thickness is very small compared with a 
representative linear dimension (L) of the body, ie. cT << L. 
In this way, the solutions obtained from the boundary layer 



Re = U x/i), (f (X) 



X 




For simplicity, assume a two-dimensional, steady 
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equations are asymptotic and apply to very large Reynolds 
numbers . 

From the Navier-Stok.es equations, 

x-direction: - 1. d P + 1. + 1. B<^xv = u d u + v ^ u 

/^ax /’ax /’ay ax ay 

y-direction: - 1_ d P + 1_ a6~vy + 1_ a<Txv = u a v + v d v 

/° dv /^ay /^ax ax ay 

and the continuity equation is: 

a u + _a_v =0 (2.3) 

ax ay 

Inserting the "typical" ( order -o £-magnitude ) values, replace 
the dependent variables as follows; 

u Uo , a u ue a u , ue 

tr Co — y Co 

ay / ax i 

Then Eq. (2.1) can be expressed as 

- i_ a p + 1. a <^xx + i_ a(Txv = u d u + v a u 
/^ax/^ax /’ay ax ay 

u/ Q~xx//’ tTxv//’ ue^ u/ (2.4) 

11 111 



where the typical values were written below the terms to 
which they correspond. 

Considering turbulent flow, all stresses are of the 
same order ie. a general stress must be of order 
/I, then Eq. (2.1) becomes 



- 1. a P + 1_ d(Txv 

/* a X z’ ay 



1 






= u a u + V a u 



(2.5) 



a X 



a y 
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where indicates a quantity of the order of magnitude 

of f/l. 



Considering a laminar flow of Newtonian viscous 



fluid. 



6" XX = d u (T xy =M f ^ u + ^ v^ (Tyy = 2j / ^ ^ 

d X, (ay 3x / ^y 

In the (Txy term, dv/dx is smaller than au/ay. 
Therefore, Eq . (2.5) becomes 



( 2 . 6 ) 



- 1 d P + u 

ax ^ a/* 


1 + 


°(ff} 


= u a u + 

3 X 


V da 

a y 


(2.7) 


Similarly, Eq . (2. 


2) can be 


written as: 






-1 ^P + 1 vv + 


1 


a(Txv 


= u a V + 


V a V 




/"ay P ay 


P 


3 X 


3 X 


a y 




1 ap / 

p ay 1/ 


Ue 


Ue</* 


\ u// 


ueV 


(2.8) 


u 




I 1^ 






If we write all the viscous 


terms together 






1 a P ugV v’ 


(>/ 


1 




ue^/ 


u// 


(2.9) 


/° a y 1-^ Ug 1 


1 


1 ' J 


1-" 


1^ 





It is known that J^/ug 1 is laminar flow, so that the 

viscous terms are also of order ug <f /I ie. dP/ay is of 
order u/ //I* , but the pressure difference between y=0 and 
y= (f is of order fiug /l^and the difference in 3P/ax will be 
negligible compared to the external stream dynamic pressure, 

X 

1/2 fUg . Therefore, for practical purposes. 



££ = 0 ( 2 . 10 ) 
dy 
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For this case, since changes in P must be of the same order 
as changes in, y, the pressure does not change significantly 
through the boundary layer . 

Thus, the entire equation of motion in the 

y-direction may be dropped from further considerations. In 

this way, the following simplified equations are left for 

the analysis of a boundary layer: 

au + Av =0 (2.3) 

d X ay 

a 

- i_ a p + a u = u a u + v a u ( 2 . 1 1 ) 

P 5 X d d X d y 

a P = 0 (2.10) 

a y 

These relations are known as Prandtl's boundary layer 
equations. Unless one encounters very high Mach numbers, 
the above orders of magnitude are not changed when 
compressibility effects are considered. 

B. LAMINAR AND TURBULENT BOUNDARY LAYER 

The low viscosity fluid flow past solid bodies should be 
considered as two regions. One is the thin region near the 
boundary in which the effects of viscosity are concentrated 
and the other is the region further away from the wall in 
which the influence of viscosity is so small that it can be 
neglected. Thus, it can be stated that all viscous effects 
are concentrated in the thin region which is known as the 
boundary layer. This boundary-layer type behavior requires 
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high Reynolds numbers. Generally, the thickness of the 
boundary layer (/) is defined as that distance from the wall 
where the velocity (u) differs by 1% from that ( Ue ) 
calculated by the ideal flow analysis. 

Consider a constant-property, steady, two-dimensional 
flow past a flat plate. If u/Ue is plotted against a 
dimensionless y-coordinate, T = (Ue/vx) / , the velocity 
profiles are geometrically similar and reduce to a single 
curve for a laminar boundary layer flow. This is well known 
as the Blasius profile. The geometrical similarity is 
maintained regardless of the Reynolds number of the flow or 
of the local skin friction. For a turbulent boundary layer 
flow, since the viscous-dependent part and the remaining 
Reynolds-stress-dependent part of the profile require 
different length scaling parameters, there is no choice of 
dimensionless y-coordinate that leads to the collapse of the 
complete set of velocity profiles into a single curve. 

The conspicuous difference in profile shape between 
laminar and turbulent shear layers can be found in the wall 
flows. Because of the constraint on eddy size in wall 
flows, the efficiency of turbulent momentum transfer 
decreases rapidly near the wall. But, the efficiency of 
viscous momentum transfer is not dependent on y distance in 
the flow which has no heat transfer. 
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C. TURBULENCE MODELS 



The unsteady continuity and Navier-Stok.es Equations are 
valid in both laminar and turbulent flows. But, it is too 
difficult to deal with the instantaneous properties in 
turbulent flow because the turbulent flow has a complex 
time-dependent behavior. Thus, the following turbulence 
models are used to make the analysis of turbulent flow more 
convenient . 

1 . Prandtl's Mixing-Length 

Consider two adjacent stream layers of fluid which 
move with different velocities. If a particle of fluid 
moves from one layer to the other, a momentum change occurs 
between two layers. The fast particles which enter the slow 
layer make it faster and the slow particles which enter the 
faster layer impose a drag on it. 

The mean velocity of a stream layer is u, and that 
of the other is u + 1 ^ u/gy where 1 is the distance between 

two layers. Also, the fluctuating velocity in the 
x-direction is u', and that in the y-direction is v'. 

Prandtl assumed that the turbulent fluctuations are due to 
the difference in the mean velocities of the two layers. So 
u' = 1 du/ay ie. the fluctuating velocity in the x-direction 
is of the order of the difference in the mean velocities of 
two layers which have a distance 1, where 1 is the mixing 
length. Prandtl also assumed that all components of 
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fluctuating velocity at a given point are of the same order 
of magnitude. Thus, v' can be defined as v' = kl au/ay 
where k is a constant. 

The turbulent shear stress due to momentum exchange 
between two layers is the rate of momentum transfer per unit 
area. Then the mean turbulent shear stress on the fluid is 
-^u'v' where u'v' ®^|l du/ay| | kl au/ay| • Since the 
values of 1 and k are unknown, combine these two unknowns. 
Then, Tf = /^ 1 pu/ay|au/ay where 1 is called the mixing 
length. 

2 . Cebeci-Smith Model 

With the eddy viscosity concept, the momentum 
equation for 2-dimensional laminar and turbulent boundary 
layers can be written as: 

(bf") ' + mjtl ff" + m[l - ( f ' f] s x(f \ 

2 ^ dx dx J 

where x is the transformed x-variable 

£ = ^ m/ y 

ak 

(1 + t) (l+£) 

Let the turbulent boundary layer be a composite 
layer consisting of inner and outer regions. Then, the 
eddy-viscosity formula for the inner region is: 



( ); = 1 ‘ 




(0 i y i y^ ) 



where 1 is the mixing length for 2-dimensional flow 
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tr 







- -J. / 


X - 1 


L 1200 ^ 


^tr ' -J 



is an intermi ttency factor for a flat plate. 
= Up Xff/i} is the transition Reynolds number. 
The eddy-viscosity formula for the outer region is: 



(^/n). - c/ 



a 



(ug - u) dy 



ftr L )l L(f 'I 



where R^ >_ 5000, the universal constant o< = 0.0168. 

R^ < 5000, cC varies with Rg according to the 
the empirical formula 



o< = 0 .0168 1 . 55 

1 + A 

A = 0.55 |_1 - exp (0.243 K - 0.298 fv)J 
fl = Rq /425 - 1 

But this model is not used in Cebeci's interactive computer 
program which will be presented in Chapter IV. 



D. TRANSITION 

The boundary layer with a finite thickness starts out as 
laminar at first in the flow past an airfoil. However, the 
boundary layer becomes unstable and all small disturbances 
begin transition to the erratically unsteady condition which 
is known as turbulence. 

In the boundary layer on blunt bodies, transition makes 
the point of separation move downstream which decreases the 
width of the wake. There is an abrupt change in the drag 
curve of a sphere. This change is due to a boundary layer 
effect and is also one of the transition phenomena. 
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The process of transition on a flat plate at zero 
incidence shows a sudden increase in the boundary layer 
thickness. Furthermore, transition involves a great change 
in the shape of the velocity distribution and a large 
decrease in the ratio of the displacement thickness to the 
momentum thickness. Also, it causes a large change in the 
skin friction. 

As described in the above, the transition from laminar 
to turbulent flow plays a very important role. Since the 
transition is not an instantaneous process and the flow is 
intermittently laminar and turbulent over a certain length 
of the airfoil, a certain point where transition takes place 
cannot be assigned. At present, there are no exact methods 
to calculate the flow within the transitional region. 
Nowadays, however, two methods - Michel's method and the e 
method - are used to predict the transition empirically. 

1 . Michel ' s Method 

Michel investigated many kinds of data for 
incompressible flows without heat transfer and found this 
empirical correlation. 



Rp = 1. 174 



1 + 22400 



R, 



t>. 4-6 



tr 



where R^ = Ue 0/P, R;, = Ue x/p . 

Since the momentum thickness grows more rapidly when 
the pressure gradient is positive, Michel's equation 
involves the effect of pressure gradient. However, the 
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effect of surface roughness, which is also very important, 
is not included. 

2 . The e^ Method 

This e method is of spatial amplification theory 
based on the solution of the Or r -Sommer feld equation which 
forms the point of departure for the stability theory of 
laminar flows. 





g' 


= a = esp ) 


' X 


(^i - 


d o<r] dX 




g' 


J 

A 






d /5r / 


where 


X = 


X//*, o(i = 


A, 

o( ; 







g'(x,y,z,t) is a typical disturbance function. 

o(. = 2 If ( A is the wavelength of the disturbance). 

A 

r is the circular frequency of the partial oscillation. 

/^r is the amplification factor which determines the 
degree of amplification 

The basic assumption is that transition begins at 
the point where a small disturbance introduced at the 



critical Reynolds number is amplified by a factor of e' 



E. LAMINAR BOUNDARY LAYER CALCULATIONS 

1 . Similar Solutions 

a. Blasius Solution For a Flat Plate 

Assume - a flat plate at zero angle of attack 

- steady, 2-D, laminar, incompressible 
flow . 
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- constant viscosity 

- negligible body force 

Since the pressure along the plate is constant, 

there is no pressure gradient. The simplified Navier-Stokes 

Equation is: 

u 3 u + V 3 u = 

3x ay 

(B,C) y=0; u=0, v=0 

y = 03 ; u = U«, 

To transform from the partial differential 

equation to an ordinary differential equation, define the 

following transformation parameter: 

”7 “ / Ua > y where % - f(x,y) 

Jl/x 






( 2 . 12 ) 



4^ ~ Jl) xUflj f where f = f(^) only 

The stream function was defined to satisfy the 
continuity equation. 



n = l± =JVxv7 All! = where f = ^ (2.13) 

ay d7 3y 



3 U 

a X 


= 


df' M 
d'7 3 X 


= 


f" y 




- 1 

2x^^ _ 








- _ 


Uo»f” ^ 












(2.14) 






2x ^ 














V = 


- 3^ 


= -Jl/xUa, 


df 


37 


- 


(1/2 X 


) f 






a X 


d 7 


3 X 










= 


1/2 J VU» 


( f ' 7 


-f ) 








(2.15) 


a u 


= UCD 


df' 3'^ 


= 


. f’/ 


Ua> 






(2.16) 


a y 


d^ 3y 


M/ / 


1/ X 
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d g = u<> / 



df" = 



i 

Ua> 



f " 



(2.17) 



^V*x 

Substitute Eqs . (2.13) ' (2.17) into Eq . (2.12), then 

f£" + 2 f" =0 

(B.C) T=0; f'=0, £=0 

^ =a> ; f ' = 1 

The solution of this Blasius equation is presented in 
"Boundary Layer Theory" by Schlichting. From the 
transformation relation. 



/ = 


5 


. 0 = 5.0 


/ 




1.7208 


X 


ct, X/y v/ Rx ^ 


X 






/ V 1 = 


0 


. 8604 


9 


= 


0.664 


1 u la, 






X 






II 


0 


.332yaU^ I 


Voa 












)7 X 






o 

II 


0 


. 664 


















b 


Falkner-Skan 


Me thod 








The Falkner-Skan transformation is for 2-D, 


axisymmetr ic 


laminar flow 


, The 


simplified Navier-Stokes 


Equation 


is : 










u 


a u 


+ V ^ u = 


- 1 


_i£ 


+ w a'^u 




d X 


d y 




d X 


a 


(B,C) 


y 


= 0 ; u = 0 


, V = 


0 





Y = 06 ; u=U(x) 



Take the same as Blasius' but different with a function 
f + f(x,"^) and follow the same procedure as Blasius' using 
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- _d_P = U ^ 

/* Jx dx, then the transformed equation is: 



f"' + m + 1 ff" + m r 


1 - (£■)*] 


= X f f ' ^ f ' 


- f " a f \ 


2 *- 


j 


a X 


a X y 



where m is a dimensionless pressure gradient parameter 

defined by: 

m = X ^ 

U dx 

(B.C) -^1 = 0; f'=0, f=0 

'I = a> ; f ' = 1 

For similarity in 2-D, laminar flow, assume f is a function 
of only. Then, 

f" + m + 1 ff" + m 1^1 - ( f ' = 0 (2.18) 

(B.C) = 0 ; f' = constant, f' = 0 

= to ; f ' =1 

The fact that m is a constant leads to: 

U = c x^ 

where c is also a constant. 

In the case of m = 0^ ie. U is a constant, 

Eq . (2.18) reduces to the Blasius equation. I f m =• 1 which 

means a 2-D stagnation flow, Eq . (2.18) becomes the Hiemenz 

equa t ion, 

fir r + f f » _ ( f ' )^ + 1 = 0 

Some solutions of the Falkner-Skan equation for various 
values of m are presented by Cebeci and Bradshaw. 

2 . Integral Methods 

a. Integral Momentum Equation 
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For steady, 2-D and incompressible flow. 



a 


^ a + 


V a u = 


1 


a p + ii a u 


(2.19) 




a X 


3 y 




ax a 






^ a 


+ ^ V = 


0 




(2.20) 




ax 


d y 








From Eq . 


(2.20), 


u / a u + 


3 V 


\ = 0 






V a X 


a y 






At the edge of the boundary 


layer , U ( x ) a u ( x ) 


= - i ap 










a X 


dx 


Then, Eq. 


(2.19) 


becomes : 








A u^ + 


^ ( uv ) 


= U(x) d U(x) 


+ yj a'^u 




A X 


a y 


dx 




ay-" 




Integrate 


this equation with respect to y, from 


y = 0 to 


y = / , using T 


= ^ au 












ay 




/ 




/‘f . - 






( 


Tw 


1 u dy 


- 


U ( X ) a u dy 


- j 


U ( X ) dU ( X ) dy 


^ (2.21) 




4 


? X 




> dx 



Also we know, / = f /l - u \ dy 

Ja ^ U ( X ) / 



U( 

0 = u f ^ ■ u \ dy 

hU(x) V U(x)y 



( 2 . 22 ) 



Substitute Eq . (2.22) into Eq . (2.21), then, 

^ = r 



dU (x)0 + J" U ( X ) dU ( X ) 

dx dx 









This equation is known as the momentum-integral equation of 
boundary layer theory, or as von KArmAn's integral equation, 
b. Pohlhausen's Method 

Assume a polynomial of the fourth degree for the 
velocity function. 
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-i 3 4- 

u = f(A) = aA+ bA + cA+ d-A 

U(x) 

where A is the dimensionless distance from the wall, 

A = y// and 0 < A < 1 

(B.C) A = 0; f = 0, 

A = 1; f = 1, f = 0, f" = 0 

Also U(x) d U ( X ) + ^ ^ u = 0 when A = 0 

dx 3 y* 

/ d U ( X ) = - 2b = 

y dx 

where J\ is the shape factor. 

then, f(A) = F(A) + ilG(A) 

3 4 

where F(A) =2A-2A + A 

G(A ) = 1/6 A (1 



Thus the boundary layer parameters if, $ and Tkv can be 
determined, if the velocity profile is known, as follows: 




/ 10 120 

Q / 37_ - A. - A: ] 

if 63 I 5 15 144 ' 

Tw = ^Ul2Ll f 2 + A 

c. Thwaites' Method 

The integral momentum equation can be written 



as : 



d© 


+ e (H 


+2) dU(x) 


= 


C£_ 




(2.23) 


dx 


U(x) 


dx 




2 






(B.C) y 


= 0 ; d*u 


= - U(x) K , 


^ u 


= 


U(x) 


L(K) 




>1 


0‘ 


ay 




9 
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where H 



^ L 
d U( X ) 



e/ = 
2 

K = 






£_ 

e 



(=> U^(x) 

X 



0 



dU(x) 

dx 



then Eq. (2.23) can be reduced again. 



U(x) d 
y dx 



= 2 



j L(K) - [h(K) + 2 ] j = F(K) (2.24) 



Thwaites writes an expression for F(K) 
F(K) = 0.45 - 6K 

Then we can write Eq . (2.24) as 

J. 



u(x)e 



0 . 45 



/ ^ J- 

I U ( X ) dx 

f ft 



U-^(x) ^0 

If Q is calculated for a given external -vel oc i ty 

distribution, H and Cj can be determined with the following 

relations. 

For 0 i K i 0 . 1 

L = 0 . 22 + 1. 57K - 1. 8K 

H = 2.61 - 3.75K - 5. 24K 

For -0.1 i K 10 

L = 0.22 + 1.402K + 0 . 018K 

0.107 + K 

0 . 0731 + 2.088 



H = 



0.14 + K 

3 . Finite-Difference Methods 

The finite difference methods are the most flexible, 
practical and efficient methods to solve the boundary layer 
equations . 
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The box method, presented by Keller and Cebeci, is 
introduced here as the preferred finite-difference method. 



The momentum equation achieved by the Falkner-Skan 
transformation can be rewritten in terms of a first-order 
system of PDF's. Then the resulting non-linear system is 
linearized by Newton's method. Finally, the block 
elimination method is used to solve the linearized 
difference equations of the boundary layer problem. 



Figure 2.1 Net Rectangle For Difference Approximation 
where | = x and = 1 + ^t/ V 



a. Box Method. 



Using the new coordinate system, | and r^, the 
transformed momentum equation for steady, 2-D, 



incompressible flow becomes: 



f = u 



(2.25a) 



u ' = V 



(2.25b) 



(bv) ' + m 4- 1 fv + m( 1-u ) 

2 




V 




(2. 25c) 




K 



n 
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A prime denotes differentiation with respect to "f . 

The boundary conditions are 

7 = 0; f(l, 0)=0, u(^,0) =0 (2.26) 

7=7 ; u ( ^ , 7^ ) =1 

We denote the net points shown in Figure 2.1 by 

,% = 0; t = i 



+ Krt n = 2, . . . , n 

fo " ° ^ j ^ ^ ' 2 , . . . , J 



(2.27) 



^7= 

Here n and j indicate sequence numbers. 

With = g( ^ n, T j ) denoting the value of any 
quantity at the mesh point ( f n, ^j), centered quantities 



can be written as two-point averages: 

^ = l/2( ^ + I ), l/2( 7j + ^j-1) 

n fl-» n n ^ n 

= l/2( f. a J; t ? ) 



( 2 . 28a) 
( 2 . 28b) 



n 

The finite difference approximations of Eqs . (2.25a) 



and (2.25b) for the mid point of the segment P, , are: 

(2. 29a) 



n 



.c /I 

f .f-t = u 






n 






\/ 



(2. 29b) 



Eq . (2.25c) can be approximated similarly by 

centering about the mid point of the rectangle P1,P2,P3,P4 

1. Centering Eq . (2.25c) about the point 

.n-’/i 

< I ' 7 > 

Let the left-hand side of Eq . (2.25c) be L and use 

Eq. (2.28b) 
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n n-i y 

1/2 (L + L ) = I 



n^'/i 



n-Yj. 



n \ n-Vj.! 0 



2K, 



i n j. n-i n n «-/ /> 

(u ) - (u ) - (fv) + f V - V f 4- (fv) 



/)-' fl n-/ /> 



n- t 



By rearranging, 

(bv')” + m + 1 (fv)'’ + m' 



2 

,n-'4 



[l - (u*)"] 



i Pa'’ 

= - L + -2 ) - (U ) 



-(fv) + f v-v f + (fv) J 






M 






rt-i J. ^ 






(bvM r ^,(fv) -cx^(u ) + <?<(v" f - f v") = R (2,30) 

where c< = ^ - . . . , . >f 



±n-!/jL ^ 

__Z . ( 7 < = m ^ 1 + • cK = m ^ a< 

K. ^ 2 






/7-I 



R =-L +c?(|(fv) -(u ) 

n -I 



[< 



r?'l 






- m 



[' 



(bv)' + m 4- 1 (fv) 4- m(l - u'*') J, 

2 



n - ) 



( I 7 • ) Using Eq . (2.29), 

7 ^ ^ 



2. Centering Eq. (2.30) about the point 






r l" n n In 

I ( bv ) ' J . = bj vj - b».)Vj/-i 

hj- 



The Eq . (2.30) becomes 

h'l (b} v/ - hi, ^1^ ) 

i»-i _ r) _ ""I n _ 



4- o( ( V. , f V- “ f ■ , V - ,. ) = Rv V 

i'A i'X. i-A 



where 

D ^ * — 

%--)i ■ 


[- cK 


(fv)^“,;£ - 


l""' = 


' 1 

rr 

•A, 


- b,-., vj:., ) 






n-i 



The boundary conditions at 



2 

n 



- m 






(2.31) 



.-U^)^--J4.j 



n-i 



, are 



■ o 
ft 

'7 



0 , 

1 , 



u;' = 0 



(2.32) 
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b. Newton's Method 



For simplicity, let ( f , u^- , 
( fj', Uy , Vy ) at ^ = I . 



) be 



Then Eq. 


(2.29) and 2.31) can be 


written as 








- h^; (uy f 

2 


uy., ) = 0 


( 2 . 33a) 




- '^3-1 


- _hi (vy + 

2 


vy_, ) = 0 


(2. 33b) 


h‘. 

3 




- by.| vy_i ) + 


ek, ( fv)y_;i - (u"^ 


h'-'/L 






-1 n-l 

-'A ^3-'A - ^ i-y. 


rt 'I 

vy.;i ) = Ry,^ 


( 2 . 33c) 


Here the 


unknowns 


n - 1 

are on the left-hand side and R y - 




involves 


only known quantities. 








Now, 


Newton's Method 


is applied to turn 




Eq. ( 2 . 33) 


into a 


linear system complemented by boundary 


conditions Eq . ( 2 . 


32) and initial 


values 




i‘:‘ = 


0 


c<7) 

u, = 0 


fl-i 

Vo = Vp” 




= 


fi-i 


(<r) n-' 

u •. = U • 

3 3 


Vj. = Vj. (1 <J < 


J-1) 


^7 = 




(») , 
1 


II 

< 




The superscripts 


in parenthesis 


represent the iteration 


number as 


follows 


I 






f 


(0 

/ ' '^3 


0'> 

, Vy 1 


=0, 1, 2, ... 





A>*') -h’) f' o> 0-*\) 0> - (I) (^,'+0 {)) ^ ((' 

^3 = + /£^- . Uj- = a^. t/Uj. , Vj. = Vj. + /Vj. 

where /f << f, /u << u, /v << v. 



Replace fj, uy, vy in Eq . (2.33) with these expressions 

Then Eq. (2.33a) becomes: 



r r . r <•> 

- ^^'-1 



U) 

'^3-' ^ ^ 



II) 1 
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Thus 





• - ^ 


1 

1 

^k. 


u^- + /u/., ^ 


1 = h' 




(2. 


34a) 






(1) 


in 


iiy 








where 


(r,)^- 


= f - 












From 


Eq . ( 2 . 


33b) , 












o'> 


+ /uy - 


0) 


;■ , - [' 


+ / Vj- + 


U) 

V. 


r <" 
/ V/., 


]- 




/u^- - . 




V,- 1 /v^-., ] 






(2. 


34b) 








<i; fi> 


(i; 










where 


V-.= 












From 


Eq. (2. 


33c) 












‘>V 


HP . 

( b / V 


il; (u p /o 

/ - b/-. ^vy., ) + ^ 


(, / ( fv)y. 


- lAj^f ( u 


)j->i 






+ ( 


rt-i p (O 

v,-.y^ ^ f^V 


;i - ^y'/x ^^vyx ^ 










1 r 


. <•> (M 


(i; U) 


o') 




.1 (n 


- ^3-'^ ~ L" 


,■ (b,- 


- '^r 




o^'i ( U 


V-v. 




^ ( 


^ n-l , h‘< 


# 0 -' in 

h'-'A ^3-A 


)] , 








here 


1/2 


[£”Vv;> * 


v|V£-:' 1 


in l!, 


+ 

1 


4/^^t 


/ ( u 


*)'■' - 




t 4. U- 


(o 

y- 1 








r u ; 


1/2 ( / f 


(»> 

>■ * 


) 








/v!'^ 


1/2 ( /v 


<*’> r (o 

) * “^v. 


) 









then (s, 
+ 



where 



)j,' / vj + (Sjl)v/v.-., + £.• + (s^)y /fv_, + (Sy)//uv 

/V. = (2.34c) 



(r^),' 


- r'*" - 

" ^3 - A. 


r , /!> 

h • (b . V . - 
1^33 


•^y-. V- ■" 


- 






rt-i (ip 1 

- fj-A V/^ ’ J 


V 


, . o’l 

= hj by + 


iJ. <=^'^ 

2 ^ 


- e-( 

2 


f 

^->5. 


(s,)^< 


- » ( »; 
= - by by.. 


+ o/| 

2 


f - 

3-' 


, ff'i 

4 




= v'.'" 

2 ^ 


+ 0<f 

2 


ff-l 
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oi 

~2 



(24)3 = 


0/ i 
2 


(1 > 




- 


0/, u^- 


II 


- 


1 (!) 



the boundary conditions are: 
/ft? = 0 /uo = 0 



V 



3 -'A 




(2.35) 



c. Block Elimination Method 

Eqs.(2.34) and (2.35) are the linearized 
difference equations of the momentum equation for external 
flows which has a block tridiagonal structure. We can solve 
these equations by means of the Block Elimination Method as 
discussed by Keller (1974). 

Let us define the three-dimensional vectors, /- 
and r', to express the system in matrix-vector form, 

<J 



■f. 



^^3 



0 i j i J 













0 




(r, 


0 

II 


0 


II 


(r^ V 








J 



1 i j <. J-1 (2.36) 
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= 


(r. 








j = J 










0 












and the 3 


X 3 


matrices 


, A,- 


B ' C ' 






1 


0 0 










1 -h^' /2 


0 


Ao = 


0 


1 0 








= 


(Sj)^- (S^)j 


(S, )^- 




0 - 


■1 - 


h, /2 








0 -1 


/ 2 



1 < j < J-1 



Aj = 



Si = 



1 -hj/2 0 

(Ss)j(Ss)j (S,)j 

0 1 0 , 

1 - h ^/2 0 

0 0 0 



j = J 



(2.37) 



1 < j < J 



= 



0 i j i J-1 



0 0 0 

0 0 0 

p 1 -hjVi /2_ 

Then Eqs. (2.34) and (2.35) can be written as 

^ = R 



(2.38) 
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where 



^ = 



A C <? 

B/ 



^3 



Bj-i Aj. I Cj., 

B 7 Aj 



/ = 



' io ' 




' ' 


J, 




r, 




K = 


O' 






0 '> 









Le t 


us factorize the matrix ^ 


to 


1 solve Eq . (2.38) 




^ = xy 




(2.39) 


where 








I 




Vo 




1 

1 

V 

X 




y. _ 


11 




■3^ 










yj-. ^7-. 








yj 



Here I is the 3x3 identity matrix. 

x^- and y^' are also 3x3 matrices. 
According to Eq.(2.39), we can find 

Yo = Ao 

Xj yj ., = B j ( j = 1, 2, . . . , J ) 
y^- = A^- - X J. c , ( j = 1,2..., J ) 



( 2 . 40a ) 
( 2 . 40b) 
(2.40c) 



and the matrix x^' has the same structure as that of the 
matrix Therefore, x^- has the elements like this. 



* 3 ' 



(X „ ) J 

( X i, ) 

0 



3 



( X 

( X ),• 



(X,J )^' 
(X^J ). 
0 
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and yean be denoted as: 

4 







(yn )a 


(y,x )i (y,3 








y = 


(ya, )^^ 


(yji (y^3 V 










L 0 


-1 -h,'.., 


/2^ 




From Eq. 


(2.40a) 










(y .1 )o 


= 1 


(yii )o 


= 0 (y,j )o = 


0 






= 0 


>1 


II 

^0 

>1 

r-4 

II 


0 




From Eq. 


(2.40b) 










(X (1 ), 


= -1 


(X ,j. ), 


= -1/2 h, (x,3 


)/ 


= 0 


(X 


= (s^) 




= (s^ ), + (X i^) 


/ (X,j) 


From Eq. 


(2.40c) 










(y„ ),' = 


1 (y 


iJ. V 


-Jl3 -(x,3 )/ (y ,3 

2 




II 


II 


(s, )^-(y 




= ( s 



Then we 


can compute the 


elements of x with Eqs. 


( 2 . 


37) 


and 


(2.40b) 


for 1 j 


L j < J, 










(x„ )^- = 


1 } 


<y^3 V-. 


+ M r Jy' (y^i V-,- (yx>. 
2 L 2 


^3- 1 


J] 




yo * 






■J, i 


'(f)^ * 


‘Xii V [<y,iV-. 4 ■ 


‘y,3 V-. 


11 


II 


-M 

2 


[ ( X „ ) j' 


'vy >j'-, * <='w<3 <y^jij'- 


.] 






( X )^> = 


1 ; 


j (Sa)/ (y^, 


)/-, + ( 54 )^' (y»j)a'., + 

[(svV(y;^V., 








yo 


r (Si)a'(y^, V_, 


(x^)j> = 


1 , ' 


f ^'(S^ )a 


- (s^)^‘ + (X )^' [(y,^ )^'_, 


- hi 




V.,] 1 


y. 1 


2 


(X^j)^- = 


(y.a V- 


, + (x^ )^- (y^ )j>_, - (s^ 


h' 






>1 

II 

o 

>1 




- (yi3)j'., (yn )a'., - 

(y^, V-. ■ V-. 


V-.] 
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V = _h» (y 



XJL 



V-. 



- 



3 ^'1 



Let y 


/ = z. 


then xz = R 




(2.41) 


Thus, 


Zo = 








(2.42) 




Z^' = r. 


- ‘ 


^ 3 " 


(1 i j 1 J) 






where 


N 

II 


(Z, )/ 

( z^ )y 
(Z,). 

" J 


(0 i j i J ) 





From Eq 


. (2 


. 42) 


for j 


= 0 


















(z, )o 


= (r, 




(Zx)o = 






(Z3 




= 


)o 




and 


for ; 


^ 1 j < 


. J 


















(z, 


V = 


V 


- 




- 


(X ,x 


V 




- (X 


13 




(Zx 


= 


( rx )j' 


- (X;,, 




- 


(Xjj. 


)^'(Zx 




- (X 




V ^ 


(Zj 


= 


( rj ) ,■ 




















From Eq 


. (2. 


'll). 
























i'j/j = 
















( 2 


.43a) 








^3 - 


^3 J V 1 


(0 


i j 


<. J- 


1) 




( 2 


.43b) 



Then the vectors can be calculated with Eq . (2.43). The 

three components of , for j = 0 , 1 , 2, . . . J-1 , are: 



/Uy = - J^+, /vv - e, 

2 ' 



/ = I [(y„)/[(Zi)/ + e,(y_,^)^-] - (y (z, 

/wj- =_I_ [(z,)y - (y,^ - (y,^)^- /v^-1 

(y„V 

where e = ( z 3 ).- - / uv^ + Ji^+i v v 



e ( y ) . ( y ) , 
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ii 

1 


hi 

2 


.[< 


Ya, V (y,a ^4 + 


(y^ 


>1 






- 


(y^, )/ (y,^ ) 


/ + ( 


Yaj )/ (y II 


for j = 


J 










/ Uj = 


( 










-A 

t) 






(Yii )j - ej 


(y« 


)cr 






(y 


,3 )o-(y w V - 




)j(y n V 


/ Wj = 




© A 


~ ( y (J ) J 




















wh 


ere 


= (Z| )J - 












® J “ ( ) J - 


(y^j 


)j /u J 


;e calcul 


at 


ions 


are stopped 


when 










1 < 


e 




where 




V ( o ) 


is the wall 


shear 


paramete: 



^ is a prescribed value. 



F. TURBULENT BOUNDARY LAYER CALCULATIONS 

Turbulent fluid motion is an irregular condition of flow 
in which the various quantities show a random variation with 
time and space. Therefore, turbulence is characterized by 
random and chaotic motion of fluid particles. 

The velocity varies randomly at any point in a turbulent 
fluid. The velocity components in a three-dimensional 



turbulent 


fluid 


are 


• 


u 


= u 


+ 


u ' 


V 


= V 


+ 


v' 


w 


= w 


+ 


w' 
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where u, v, w represent the mean velocities 



u = 



1 ( u dt, 

t, - 1, X 



and so on. 



u', v', w’ represent the amount of fluctuation of the 
instantaneous velocities from the mean velocities. 

Obviously, we can see that 

dt = u' = o, etc. 



1 (u' 

tx - t, 



Similarly, the following quantities can be defined, 

dt = u' , etc. 



1 ( u' 

ta - t, \ 

1 u'v' 

t, - t. 4. 



dt = u'v', etc. 



Even though the flow is turbulent, the time-dependent 
momentum equations are valid. However, it is impossible to 
solve the momentum equations for turbulent flow because the 
fluctuations are random and chaotic. Thus, Reynolds 

modified the momentum equations by introducing the mean 

U. 

value and fluctating values of the flow quantities and by 
assuming the fluctuations to be continuous functions of time 
and space. From the momentum equations. 



<»u + u_d_u + V ^u + w i^U = »^P + j) / o''‘u + ^^u +_£ji \ 

^t ay a-z 9 X (ax'* ay‘ <)zv 

By multiplying the continuity equations with u 

u I ^ u + ^ V + a w ^ = 0 
\ a X ay a z / 
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Sum these two equations and rearrange, then 



3u = 


- 1 ap 


+ J.I 


v^au - u 




^ 1- 


)>au - uv j + a 


/Pau 


~ uw 


■5 t 


1° dx 


ax \ 


ax 


I 


ay 1 


ay / azi 


L a z 




Here 


take the mean 


value 


of 


each 


term using u' = 


= 0, 






a( u + 


u ' ) 


= ^ u + 


au 


f j- 


a u 








a t 




at 


at 




<a t 








- 1 


a(P + 


P' ) = 


-1 


ap 










r 


a X 






a X 










^ a(u 


+ u' ) 


- (u 


+ 


u' 


1) 

1 


2 uu' 


/ ^ 
+ u) 




V 


ax 








ax 


















= v>aQ - u*" - 


u 








^(u_j_u_L) = 9^u 



d y 



^ u + u ' ) 
^ z 



dy 

= 

d z 



uv 



- uw 



- u'v' 



u ' w 



(2.45) 



Substitute Eq . (2.45) into Eq . (2.44) and rearrange. Then, 






ju + u <)u + 

at ax 

= -_aj_+ 

^ X 



V <)u + w au \ 

ay a z / 

yU <]^U~ f ^ *+ -5 u ' V * 



+ a u ' w 

a z 



( 2 . 46a) 



By similar procedures for the other directions. 



a ( + u ay. + V av + w ay 

' at ax ay az / 





= - a P 

ay 


+ X (?*v - al av'u' + 
"\ ax 


d v' + a V ' w ' 'j 

ay d z } 


(2.46b) 


r 


a w + u 

a t 


a w + V ^ 

ax a 


w + w 

y 


a w \ 
a z / 






= - a p 


-^1 


/ w ' u ' 


+ ^w'v' + ^w' 'j 


( 2 . 46c) 




^ Z 


1 


ay a z / 




and 


from the 


continuity 


equation. 






^ u 
a X 


+ a V + 

a y 


d W = 

a z 


0 


(2.47) 
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In laminar flow, we have 4 unknowns and 4 equations. 



But in turbulent flow, we have 10 unknowns and only 4 
equations as shown at Eqs.(2.46), (2.47). This is the 

reason why we need the turbulence models which were 
discussed in Section C. 

G. SEPARATION 

In some cases, the boundary layer thickness increases 
considerably in the downstream direction and the flow in the 
boundary layer reverses its direction. The change in 
direction causes the decelerated fluid particles to move 
outwards, which means that the boundary layer is separated 
from the wall. This phenomenon, boundary layer separation, 
is always related to the formation of vortices and to large 
energy losses in the wake of the body. 

Let us consider the simplified boundary layer equation 
in order to investigate two-dimensional separation; 



u j u 
dx 



+ V u 

dy 



1 P + u 



ax 



and define the separation as 




0 



Since u=v=0aty=0. 
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The velocity profile in the boundary layer always has an 
inflection point in the region of decelerated flow and the 
velocity profile at separation point must have an inflection 
point. Thus, separation can occur only when the flow is 
decelerated. In two-dimensional separation, a bubble of 
fluid which has low velocity is always formed. This bubble 
is often unsteady and distinguished by interior streamlines 
which are closed loops or which extend to infinity. 

There are two severe problems in boundary layer 
calculations when separation occurs. 

1. The Goldstein singularity at the separation point in 
direct boundary layer calculations. 

2. Numerical problems downstream of the separation point. 

If a boundary condition prescribes the pressure 

gradient, then the boundary layer methods suffer from a 
singularity due to separation. Moreover, the singularity 
causes the numerical breakdown in direct boundary layer 
methods near the separation point. This Goldstein 
singularity can be overcome by using the displacement 
thickness or the wall shear stress, instead of the pressure 
gradient, for the boundary condition, which makes it 
possible to integrate the boundary layer equations through 
the separation point. Also, the full Navier-Stokes 
equations exhibit no singular behavior. 
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On the other hand, the boundary layer equations lead to 
a numerical instability in the region of reversed flow. The 
FLARE approximation is the most common method to overcome 
this instability. The momentum transport term 
u ^u/<?x is deleted where u is smaller than zero. But the 
accuracy of this approximation decreases as the region of 
reversed flow increases. Therefore, it is necessary to 
introduce an upstream influence in that region. This is 
incorporated by the downstream upstream iteration procedure 
(DUIT) which consists of a sequence of alternating up- and 
downstream sweeps with the momentum transport term u <^u/dx. 
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Ill . 



INTERACTION METHOD 



A. INTRODUCTION 

If the flow remains attached and the Reynolds number is 
high, the pressure distribution and the overall lift force 
can be obtained from a potential flow solution. 

Conventional boundary layer methods provide additional 
information about the skin friction distribution and the 
overall drag force. However, if the flow separates, no 
information is available for regions downstream of the 
separation point. Therefore, the overall forces cannot be 
obtained . 

For this reason, interaction methods are introduced to 
overcome this problem. They provide a special coupling 
between the inner viscous and the outer inviscid flows, and 
can be classified as follows: 

1. Weak interaction methods 

a. Direct method 

b. Inverse method 

c. Semi-inverse method 

2. Strong interaction method 

The weak interactions provide only a loose coupling 
between viscous and inviscid regions, i.e. two different 
regions are treated alternately. The viscous flow solver 
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deals with the flow in the viscous region and yields the 
boundary condition of the inviscid region, while the 
inviscid flow solver deals with the flow in the inviscid 
region and yields the boundary condition of the viscous 
region. The exchange of information through the boundary 
condition is slow, but regions such as those near separation 
and trailing edges require fast and direct coupling between 
viscous and inviscid region. While the strong interaction 
method treats displacement thickness and external velocity 
simultaneously, the weak interaction methods process one of 
these quantities as input and the other as output. 

1. Inviscid flow methods 

a. Direct boundary conditions 

(1) Prescription of the airfoil shape 

(2) Zero normal velocity at the surface 

b. Inverse boundary conditions 
Prescription of a velocity distribution for 
the unknown airfoil shape 

2. Viscous flow methods 

a. Direct boundary conditions 

(1) No slip condition requiring zero normal 
and zero tangential velocity at the 
sur face . 
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(2) Prescription of the external velocity 
i.e. u-component of velocity at the 
edge of boundary layer 

b. Inverse boundary conditions 

(1) No slip condition 

(2) Prescription of the displacement 
thickness 

c. Boundary conditions for simultaneous 
interaction 

(1) No slip condition 

(2) Prescription of a linear combination of 
displacement thickness and external 
velocity 

B. WEAK INTERACTION METHODS 
1 . Direct Method 

A direct inviscid flow solver is combined with a 
direct viscous flow solver (see Figure 3.2). The boundary 
conditions are: 

u (x,o) = 0 v(x,o) = 0 (3.2) 

u ( X, y ) = Ue ( X) 

This method is terminated at the point of vanishing 
skin friction. In the direct scheme, the pressure is 
calculated from the inviscid region, but the displacement 
thickness is determined from the viscous region. Because of 
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this phenomenon and Goldstein singularity, this method is 
not appropriate for flows with strong interference effects 
between viscous and inviscid regions. 

2 . Inverse Method 

This method consists of an inverse inviscid and an 
inverse viscous flow solver (see Figure 3.3). The boundary 
conditions of the inverse boundary method are: 

u (x,o) = 0 , v(x,o) = 0 (3.2) 

= u(x, y^ ) [yg - 

where ^ is the stream function. 

The roles of displacement thickness and external 
velocity distribution are exchanged in this method. The 
troubles related to the Goldstein singularity can be 
overcome, but the whole procedure takes a long time due to 
very slow convergence. Thus, this method can be applied to 
the regions of separated flow only and needs severe 
under -relaxation . 

3 . Semi-inverse Method 

A direct inviscid flow solver is combined with an 
inverse viscous flow solver (see Figure 3.4). The input is 
the displacement thickness and the output is the external 
velocity distribution for both solvers. The two external 
velocity distributions are combined and then an updated 
displacement thickness is obtained through a relaxation 
procedure. A formula for satisfactory convergence is: 
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^ old ( X ) 



(3.3) 



d new 



(X) 



r 1 + w / 


Uev(^^ -1) 


L 1 


Uel(x; /_ 



where w is the relaxation parameter. 

The numerical weaknesses can be overcome, but the 
coupling is still loose. 



C. STRONG INTERACTION METHOD 

The simultaneous method solves the boundary layer 
equations subject to an interaction law (see Figure 3.5). 
Viscous displacement effects are allowed to cause 
substantial changes in the external velocity distribution. 
Since both displacement thickness and external velocity are 
treated as unknowns, one more additional relation, the so 
called interaction law, is needed. Thus the boundary 
conditions are: 

u(x,o) = 0 v(x,o) = 0 

4^ ) = u ( X , yg ) [ “ /^( X )] (3.4) 

u (x,y-) = u/(x) + l/iT f d fu (I , v> ) /tl)1 di 

I dj^ -J X -i 

The last equation represents the interaction law. 

Here, the external velocity for the boundary conditions 
can be written as: 

U (x,ye) = Ug (X) = Ug (X) + AUe (x) (3.5) 

o 

where Ug (x) is due to inviscid flow past the airfoil 
A\}ff (x) is the perturbation velocity due to the 
displacement effect of a boundary layer. 
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To obtain 4Ug. ( x ) , the blowing velocity concept is used. Let 
us consider an airfoil along which sources are distributed 
(see Figure 3.1). This surface distribution of sources 
represents the effect of boundary layers. There are two 
ways to predict the displacement effect of a boundary layer 
on the outer inviscid flow: 

1. Compute the inviscid flow past a displacement 
body 

2. Replace the condition of zero normal velocity on the 
surface with a condition of prescribed blowing 
velocity at the surface 




The streamlines are displaced away from the surface by 
the distributed sources which eject the fluid at the 
surface. Thep the virtual displacement body becomes a 
streamline and the flow tangency condition is: 

V (X, / ) = J</ (3.6) 

Utf(x) Jx 

From the thin airfoil approximation, the displacement 
thickness is assumed to be so small that u-component of 
velocity do not change across the layer and the airfoil in 
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this connection can be represented by a straight line. 
Therefore, the blowing velocity v(x,o) is equal to half of 
the source strength <T(x). 



V (x,o) = 1/2 (T (x) 



(x,o) = v(x,/ ) 



= d/ 



dx 



1 ^ V dy 

f^dUe 

dx 



/*) 



= d (U, 

dx 

Thus, d (.UeS^) = 1/2 <T ( x ) 

dx 



(3.7) 



(3.8) 

(3.9) 



Now, the perturbation velocity due to the displacement 
effect can be written as: 

/'Ufi(x) = 1_ [ ^ Slh (3.10) 

27 T j X -i 

'' 1 

Finally, Eq. (3.6) becomes: 

Ug(x) = Ue (X) + ^ ( (T( >) d> (3.11) 

27 T j X - ^ 

*1 

Eq . (3.11) represents the interaction law in usable form and 

the integral on the right hand side is known as the Hilbert 
integral . 

In the simultaneous method, the inviscid flow solver 
provides an initial external velocity distribution, but the 
inviscid flow solver is not incorporated in the overall 
iteration process. Thus, the viscous flow solver needs the 
interaction law to obtain the necessary information about 
the inviscid region. 
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Figure 3.4 Semi-Inverse Method 
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Figure 3.5 Simultaneous Method 
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IV. DESCRIPTION OF CEBECI * S INTERACTIVE PROGRAM 



A. INPUT DESCRIPTION 
1 . Introduction 

The Reynolds number is the most powerful parameter 
which can affect the flow. Wholly different flows can 
result from different Reynolds numbers (e.g. 6.0E + 06 
against 0.28E + 06). Even though a flow with high Reynolds 
number behaves well and does not separate, it may have an 
extensive flow separation at a low Reynolds number. 

Numerical breakdown in boundary layer computations is partly 
due to this extensive flow separation at low Reynolds 
numbers. To avoid unrealistically large regions of flow 
separation, it is necessary to make some changes in the 
turbulence model for transitional flow. Reduction in 
transitional region to a proper level is a method to 
decrease numerical problems in computation. By doing this, 
more reasonable results can be obtained in low Reynolds 
number flows. Also, the computation at low Reynolds numbers 
needs more iterations to converge than at high Reynolds 
number, and the lift coefficient increases with Reynolds 
number at the same angle of attack because low Reynolds 
number flows exhibit a stronger displacement effect than 
high Reynolds number flows. 
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The results may be affected by transition location. 
The transition location can be either fixed by the user or 



computed from the following empirical formula given in 
Cebeci and Bradshaw (1977). 



Rg = 1. 174 



M + 22400 \ 

^ R, ) 



p. 4 6 



(4.1) 



Laminar flow separation, however, may occur upstream 
of the transition location predicted by the above equation. 
In this case, it is assumed that the onset of transition 
corresponds to the point of laminar separation, and the 
following message will be issued in the output: 



TRANSITION LOCATION HAS BEEN RECOMPUTED AT THE 
POINT OF LAMINAR SEPARATION 



The results of this program agree well with the 

experimental results up to stall at high Reynolds number 

flows. However, low Reynolds number flows may not agree 

well or experience numerical breakdown close to stall with 

the following messages: 

at the very top of output, 

+ IFY 002 I STOP 1 

or at the very bottom of output, 

IFY 207 I VFNTH: PROGRAM INTERRUPT (Z)- 

FLOATING-POINT EXCEPTION OVERFLOW 

IFY 259 I STNCT: /ARG/ = /argument/, 

(HEX = hexadecimal), APPROACHES SINGULARITY 
or MAXIMUM NUMBER OF ITERATIONS EXCEEDED 

Thus, the range of angle of attack should be considered 

carefully . 
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If the user requires many sweeps with IPRNT = 0 or 
2 , then the user needs to choose class "G" or "J" for enough 
running time and add additional command ", LINES = (m)" just 
after nnnnP on the following line to make enough space for 
print: 



11* MAIN ORG = NPGVMl nnnnP , LINES = (m) 

where nnnn is user's ID number. 

m is the number of output lines in thousand 
(Generally 10 }s enough) 

If the user does not do that, the following error 
message will be issued at the very top of the output and 
printing will be stopped abnormally: 

lAT 1600 JOB number (user's job number) LINES EXCEEDED 
lEF 4501 user's job number - ABEND S722 UOOOO 

2 . Detailed Input Data Description 

The input to the computer program consists of 1 

title line, 3 control lines and airfoil coordinate. User 

must follow the data format for specified column and type. 

I Title Line I FORMAT (18 A4) 

This line provides any description as desired with 
any acceptable machine characters. 

I Control Line 1 I FORMAT (5 15) 
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5 

10 

15 

20 

25 



Name 



Explanation 



ITF( 1) 


Transition flag for the 
lower surface 


ITF(2) 


Transition flag for the 
upper surface 
= 1 Point of transition 
has to be specified 
by user. This option 
shall be used if ex- 
perimental results 
are available or to 
avoid oscillations of 
the computed transition 
points. In the latter 
case, fix transition at 
the most upstream 
occurrence of computed 
trans i tion . 

= 4 Point of transition 
will be calculated 
according to Eq . (4.1). 

No input is necessary 
for XTRL and XTRU. The 
computed point of 
transition will be 
redefined, if laminar 


IRSTRT 


separation takes place 
upstream of the 
location predicted by 
Eq . (4.1) 

Boundary layer restart 
flag . 


Read Starting 


Solution Store Final 

solution or Unit 2 


0 No 

1 Yes 


Yes 

Yes 



Leave 0 for general use 



IGLMAX 


Number of sweeps. 

Low Reynolds number 
flows need more sweeps 
to converge. 


IPRNT 


Print flag to control 
output print. 
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= -1 Summary print for the 
last two sweeps 

= 0 Summary print for each 
sweep . 

= 2 Whole print for all 
sweeps . 



1 Control Line 2 


1 FORMAT 


(4 ElO.O) 


Column 


Name 


Explanation 


1-10 


RL 


Chord Reynolds number 


11-20 


XTRL 


Fixed transition location 
for lower surface 


21 -30 


XTRU 


Fixed transition location 
for upper surface when the 
stagnation point is above 
the leading edge, XTRL may 
be positive. When the 
stagnation point is beneath 
the leading edge, XTRU may 
be negative. 


31 - 40 


ALPO 


Angle of attack in degree 


1 Control Line 3 


1 


FORMAT (15) 


Column 


Name 


Explanation 


1-5 


MPTS 


Number of airfoil 
coordinates . 



I Coordinate Data Line I FORMAT ( 2F 10.0) 



Input dimensionless airfoil coordinates as two columns in 
one row ( x/c, y/c) . The order is trailing edge — slower 

surface — > leading edge > upper surface ►trailing 

edge. Thus, trailing edge is input twice. 
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B. OUTPUT DESCRIPTION 



1 . Introduction 

To take the appropriate results from the output, 
user must check the convergence. This can be done by 
comparing the convergence indicators, lift coefficient and 
displacement thickness at trailing edge, when the 
computation is completed successfully for the given sweeps. 
If the convergence indicators show steady values over some 
sweeps i.e. each difference of the convergence indicators is 
less than 1%, the results are considered as converged. 

Sometimes, the user may experience failure in 
computations at low Reynolds numbers and high angles of 
attack due to numerical breakdowns. These breakdowns take 
place when Newton iteration does not converge near to the 
trailing edge on the upper surface or the computation is 
terminated by Fortran error with this message: 

IFY 2511 SSQRT; ARC = argument, LESS THAN ZERO 

Also, extensive flow separation cause numerical 
breakdowns of the boundary layer computation. These 
unrealistically large regions of separation at low Reynolds 
numbers can be reduced by changing the transitional flow 
model . 

2 . Detailed Output Data Description 

The output of the computer program for IPRNT = -1 
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and 0 can be divided into three parts as follows; 

1. Input data and inviscid lift coefficient 

2. Alternating boundary layer parameters for each 
surface and wake. 

3. Inviscid and viscous pressure distributions 



The 


following is 


the list 


of 


computed boundary layer 


parameters 


in the order 


in which 


they 


appear in the output 


print. 










Name 


Name 








in print 


in proqram 


De f ini t ion 


Explanation 


X(NX) 


X(I) 






Dimensionless sur- 






c 




face distance from 
the stagnation point 


X/C(NX) 


XC(I) 


X 




Dimensionless chord- 






C 




wise distance from 
leading edge 


V( 1,NX) 


VWL(I) 


/ d ^ f 




Dimensionless wall 






V 


fw 


shear stress para- 
me ter . 


CF 


CF(I) 


2 




Local skin friction 










coe f f icient 


DELST 


DLS(I) 


J* 




Dimensionless dis- 






c 




placement thickness. 






/ 00 




u ] dy 










Ufi / 


THETA 


THETA { I ) 






Dimensionless 






C 

(CO 




momentum thickness 




e = 




u 


'I dy 






J Ue U 

o 




1 


UE(NX) 


UE(I) 


Uel 




Dimensionless exter- 










nal velocity 



computed from the 
inviscid method with 
viscous effect. 
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W(NP,NX) 



WNP(I) 



UeV 



Dimensionless exter- 
nal velocity 
computed from the 
interactive boundary 
layer method. 



D(NX) 

DB(NX) 



IT 



D(KX) 

DB(KX) 



^ Product of velocity 

Ue V 0 \jRL and displacement 

U® c thickness in dimen- 

sionless form for 
the current sweep 
(D), from the pre- 
vious sweep ( DB ) . 

D and DB are printed 
in the output as up- 
dated values through 
the relaxation 
procedure 



f f\ey^ (N X) - 



W Vx) 
' UE{f^x) 




ITN(I) Iteration count for 

the convergence of 
solutions at a given 
NX-station . 



NP 



NNP( I ) 



Number of points 
across the boundary 
layer . 



C. HOW TO CHANGE THE ORIGINAL PROGRAM 

Cebeci's computer program is written entirely in FORTRAN 
and consists of the following eight FORTRAN files: 



FILE 


1 


FORTRAN 


A1 


FILE 


2 


FORTRAN 


A1 


FILE 


3 


FORTRAN 


A1 


FILE 


4 


FORTRAN 


A1 


T7T T TT 


C 


T7/^nT*ri ii_ 


a V 


FILE 


O 

6 


FORTRAN 


ft T 

A1 
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FILE 7 


FORTRAN A1 


FILE 8 


FORTRAN A1 


FILE 9 


FORTRAN A1 



Also, the user needs the following six JCL files to run 
changed program; 



QUEST 1 


JCL A1 


ALLPDS 


JCL A1 


QUEST 3 


JCL A1 


STOPDS 


JCL A1 


LOADMOD 


JCL A1 


INTAPCLG 


JCL A1 


"QUEST 1" 


lists all data set on MSS (Mass Storage 
System) that belong to student user ID 
number and in group PUB4B (public disk 
volume), which deletes the data set 180 
days after creation. 


"ALLPDS" 


allocates space for a PDS (partitioned 
Data Set) on PUB4B prior to reading a 
data set into it, using the utility 
program lEFBR 14 which pre-allocates or 
deletes data sets. 


"QUEST 3" 


lists the members of a PDS and calculates 
the remaining available space within the 
data set. 


"STOPDS" 


places a FORTRAN source file in a PDS. 
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The SYSUT2 DD statement describes the out- 



put data set and the SYSUTl DD statement 
describes the input data set. 

"LOADMOD" creates a load module. 

"INTAPCLG" compile^ source, code, loads and executes -text f ile. 

A list of all JCL files is given at the end of this chapter 
and more information is in the User's Guide to MVS at NPS 
( 1986) . 



The following is the procedure to change the program 
and run it: 

Step 1 

1. Make changes in FORTRAN files as required. Remember 
those eight FORTRAN files compose one program. 

Thus, the user must check all files which are 
related to changes. 

2. Review and update FORTRAN errors, if necessary. 

Step 2 

1. Change job name "xxxxxx", user ID number "nnnn" and 
library name "yyy" with user's own in all JCL files 
to make them be the user's. 

2. Submit four JCL files: 

"QUESTl" "ALLPDS" "QUEST 3" "LOADMOD". 

Step 3 



Do the following for all FORTRAN files to submit 
them to user's library. 

1. Open "STOPDS" JCL file. 

2. Change the FORTRAN file number indicated by n on the 
sixth line. 

//SYSUT2 DD DISP YYYLIB (INTAIR n) 
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3. Insert the FORTRAN file between 

"//SYSUTl DD, and 

4. Submit "STOPDS" JCL file. 

5. Check the result by any error message at the very 
top and bottom part and condition code. 

6. Update errors, if necessary. 

Step 4 

Do the following to run the changed program. 

1. Open "INTAPCLG" JCL file. 

2. Insert the input data between 

"//GO. SYSIN DD ^ " and "/*". 

3. Submit "INTAPCLG" JCL file 

From the second change, do steps 1, 3 and 4. 

D. APPLICATION OF CEBECI ' S PROGRAM 

Cebeci's interactive program was applied to a single 
airfoil, FX 63 -137, at three Reynolds numbers. The print 
flag is 0 and the number of airfoil coordinates is 49 for 
all cases. The results are compared with experimental data 
and the turbulence model has been changed to get better 
results. The experimental results are taken from Reference 
9. 

Figure 4.2 and 4.3 show the comparison of lift and drag 
coefficients for three Reynolds numbers. Drag coefficients 
are computed at the trailing edge. It is seen that the 
results are closer to the experimental data as the Reynolds 
number increases . 
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Figures 4.4, 4.5 and 4.6 show the skin friction 
coefficients on the upper and lower surfaces as a function 
of angle of attack for the same three Reynolds numbers. 
Legends indicate the angle of attack. As the Reynolds 
number increases, the length of the separation bubble 
decreases on both surfaces for the same angle of attack. 
Especially, Figure 4.4, which has a different label of 
y-axis from the other figures, indicates an unrealistic flow 
on the upper surface at low Reynolds number. In Figure 4.5 
and 4.6, however, the flow on the upper surface is separated 
from about 80% chord at all angles of attack. 

Figures 4.7, 4.8 and 4.9 present how the displacement 
thickness varies according to the angle of attack and the 
Reynolds number. As the Reynolds number increases, the 
displacement thickness decreases on both surfaces. As the 
angle of attack increases, it also increases on the upper 
surface, but decreases on the lower surface. 

Table 4.1 represents the computed and fixed transition 
locations which are used in Figures 4.10 4.13. The 

computed transition locations XTRL = 0.288922 at ALPD = 8 
degrees and XTRL = 0.206866 at ALPD = 10 degrees were 
printed incorrectly. In this case, the correct transition 
locations can be obtained by the following procedure: 

1. Calculate the difference of XTRL's at ALPD = 4 
degrees and 6 degrees. 

2. Add the difference to XTRL of ALPD = 6 degrees and 
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take the result as a temporary XTRL for ALPD = 

8 degrees. 

3. Fix transition location with the temporary XTRL 

4. Run the program to print XC(I), GAMTR(I) from 
SUBROUTINE OUTPUT in FILE 3 FORTRAN Al. 

5. Take XC(I), where the value of GAMTR(I) is finally 
equal to zero, as the correct transition location. 

6. From the same procedures with XTRL's at ALPD = 6 
degrees and 8 degrees, obtain the correct transition 
location for ALPD = 10 degrees. 

Figures 4. 10 4. 13 show the effect of variations in the 

empirical constant, ' iri turbulence model at R^ = 0.28 

6 

X 10 . In the legend of Figures 4.10 and 4.11, the three 
numbers are the values of the empirical constant, C means 
that the transition location is computed by Eq . (4.1) and F 

means that it is fixed. The lift curve is closer to the 
experimental data as the empirical constant decreases, but 
the drag curve is closest when the empirical constant is 
120 . 

The reason why we can obtain better results by reducing 
the empirical constant is shown in Figure 4.12 and 4.13. 
GAMMATR ( is the intermi ttency factor which is a 

function of the x-coordinate with values 0.0 at the 
beginning point of transition and 1.0 at the ending point of 
transition. The intermi ttency factor makes it possible to 
avoid a sudden transition from laminar to turbulent by 
smoothing out the step-shaped change of viscosity from 
kinematic to eddy. The relation between the empirical 
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constant and the intermittency factor is given as; 

3 /X 

ix - Xtr) Ug 

n ^ tr 



= 1 - exp 



i: 



dx 



Ue 



tr 

where, x^^ denotes the beginning point of transition. 

The transition length is mainly determined by the 
empirical constant which is set as 1200 in Cebeci's 

original program. Decreasing the value of the empirical 
constant reduces the transition length. In this way, the 
unrealistic flow shown in Figure 4.4 can be avoided and more 
reasonable results, which are closer to the experimental 
data, can be obtained. In Figure 4.12 and 4.13, the legend 
indicates the angle of attack, and each of four line 

patterns is used twice for two angles of attack (e.g. ; 

- 4 degrees and 4 degrees). The direction of curves, as the 
angle of attack increases, is from right to left at the 
upper surface (Figure 4.12) and from left to right at the 
lower surface (Figure 4.13). 
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List of JCL Files 



ALLPDS JCL A1 



//XXXXXX JOB (NNNN, 9999), 'ALLOCATE PDS',CLASS=A 
//)^MAIN 0RG = NPGVM1 .NNNNP 
// EXEC PGM=IEFBR1A 

//SYSPRINT DD SYS0UT=A 

//DDl DD UNIT=5330V,MSVGP=PUB^B,DISP=(NEW,CATLG, DELETE), 
// SPACE=(CYL, (^,^^,6) ) ,DSN = MSS.SNNNN.YYYLIB 

// 



QUESTl JCL A1 y^y^y^y^y^y^y^ 



//XXXXXX JOB (NNNN, 9999), 'QUESTION', CLASS=A 

//XMAIN 0RG=NPGVM1 .NNNNP 

// EXEC PGM=IDCAMS 

//SYSPRINT DD SYS0UT=A 

//SYSIN DD X 

LISTDSET GROUPCPUB'^tB) L EVEL ( MSS . SNNNN ) 

/yi 

// 



yiyiy^yiyiy^y^ QUESTS JCL A1 y^y^y^y^yiy^y^ 



//XXXXXX JOB (NNNN, 9999), 'LIST' ,CLASS=A 
//XMAIN 0RG=NPGVM1 .NNNNP 
// EXEC PGM=IEHLIST 
//SYSPRINT DD SYS0UT=A 

//DDl DD UNIT=3330V,VOL=SER=MS0005,DISP=SHR 
//SYSIN DD 

LISTVTOC FORMAT,VOL=3330V=MS0005, DSNAME=MSS.SNNNN.YYYLIB 
LISTPDS VOL=3330V=MS0005, DSNAME=MSS.SNNNN.YYYLIB 
/y(^ 

// 



y^y^y^y^yiy^yi STOPDS JCL A1 y^y^y^y^y^y^y^: 



//XXXXXX JOB (NNNN, 9999), 'PLACE PDS' 

//)^MAIN ORG = NPGVMl .NNNNP 
// EXEC PGM=IEBGENER 

//SYSPRINT DD SYSOUT=A 
//SYSIN DD DUMMY 

//SYSUT2 DD DISP=(OLD,KEEP),DSNAME=MSS.SNNNN.YYYLIB(INTAIR_), 
// DCB=(RECFM=FB, LRECL=80, BLKSIZE=8000) 

//SYSUTl DD 

/X 



// 
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LOADMOD JCL A1 






//XXXXX JOB (NNNN, 9999), ’CREATE L OADMODUL E ’ , CL ASS=C 
//XMAIN ORG=NPGVMl .NMMNP 

// EXEC FORT VCL, PARM. FORT=’LVL( 77), NOS, MOX, NOMAP’ 

//FORT.SYSPRINT DD DUMMY 

//FORT.SYSIN DD DISP=SHR, DSN=MSS . SNNNN . YYYLIBC INTAIRl ) 

// DD DISP=SHR, DSN=MSS .SNNNN. YYYLIBC I NTA I R2) 

// DD DISP=SHR,DSN=MSS .SNNNN. YYYLIBC INT A IRS) 

// DD DISP = SHR, DSN = MSS.SNNNN.YYYLIBCIfJTAIR^) 

// DD DISP=SHR, DSN=MSS.SNNNN.YYYLIBCINTAIR6) 

// DD DISP=SHR, DSN=MSS.SNNNN.YYYLIBCINTAIR7) 

// DD DISP = SHR, DSN = MSS. SNNNN. YYYLIBC I NT AIRS) 

// DD DISP=SHR, DSN=MSS. SNNNN. YYYLIBC I NTA I R9) 

//LKED.SYSLMOD DD DISP=SHR, DSN=MSS . SNNNN . YYYL I BC INTCAS ) 
//LKED.SYSIN DD DSN=NULLFILE 
// 



INTAPCLG JCL A1 



//XXXXXX JOB CNNNN, 9999), ’COMP LOAD AND GO’ 
//XMAIN ORG = NPGVMl .NNfiNP 

// EXEC FORTVCLG,PARM.FORT=’LVLC77) ,NOS 

// REGION. GO=1280K 

//FORT.SYSPRINT DD DUMMY 



//FORT.SYSIN 

// 

// 

// 

// 

// 

// 

// 



DD DISP=SHR,DSN=MSS. SNNNN. YYYL 
DD DISP = SHR,DSN = MSS.SNNJJN.YYYL 
DD DISP=SHR, DSN=MSS. SNNNN. YYYL 
DD DISP=SHR,DSN=MSS. SNNNN. YYYL 
DD DISP=SHR,DSN=MSS. SNNNN. YYYL 
DD DISP=SHR, DSN=MSS. SNNNN. YYYL 
DD DISP=SHR, DSN=MSS. SNNNN. YYYL 
DD DISP=SHR, DSN=MSS. SNNNN. YYYL 



,CLASS=G 
,NOX, NOMAP’, 



IBCINTAIRl) 

IBCINTAIR2) 

IBCINTAIR3) 

IBCINTAIRL) 

IBCINTAIR6) 

IBCINTAIR7) 

IBCINTAIR8) 

IBCINTAIR9) 



//LKED.SYSPRINT DD DUMMY 
//LKED.SYSIN DD DSN=NULLFILE 
//GO.FTOlFOOl DD UNI T=SYSDA , SPACE= C CYL , C 1 , 1 ) ) 
//GO.FT02F001 DD UN I T = SYSDA , SP ACE= C CYL , C 1 , 1 ) ) 
//GO.FT03F001 DD UNI T=SYSDA , SPACE= C CYL , C 1 , 1 ) ) 
//GO.FTOLFOOl DD UNI T=SYSDA, SPACE= C CYL , C 1 , 1 ) ) 
//GO.FT99F001 DD UNIT=SYSDA, SPACE= CCYL , C 1 , 1 ) ) 
//GO.FT08F001 DD UNIT=SYSDA, SPACE=CCYL , C 1 , 1 ) ) 
//GO.FT09F001 DD UNI T=SYSDA , SPACE= C CYL , C 1 , 1 ) ) 
//GO.FTlOFOOl DD UNIT=SYSDA, SPACE=CCYL , C 1 , 1 ) ) 
//GO.FTllFOOl DD UNI T=SYSDA, SPACE= CCYL , C 1 , 1 ) ) 
//GO. FT06F001 DD SYSOUT=A 
//GO.SYSIN DD X 

/X 



// 
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Figure 4.1 FX 63 - 137 
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Interaction method 
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Interaction method 
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Of Local Skin Friction 
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Of Local Skin Friction 
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Figure 4.7 Distribution Of Displacement Thickness 
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Fiqure 4.8 Distribution Of Displacement Thickness 
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Distribution Of Displacement Thickness 



TABLE 4.1 TRANSITION LOCATION 



ALPD 


COMPUTED 


FIXED 


XTRL 


XTRU 


XTRL 


XTRU 


-4 


0.010680 


0.584178 


0.010680 


0.584178 


-2 


0.326179 


0.542882 


0.326179 


0.542882 


0 


0.391919 


0.501232 


0.391919 


0.501232 


2 


0.458252 


0.478890 


0.458252 


0.459243 


4 


0.524927 


0.416490 


0.524927 


0.416490 


6 


0.652658 


0.353115 


0.632385 


0.353115 


8 


0.288922 


0.288922 


0.676143 


0.288922 


10 


0.206866 


0.224747 


0.718631 


0.206866 
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Figure 4.10 Comparison Of Lift Curves 
For Different Empirical 
Constants ( Re = 0.28 x 10 ) 
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Figure 4.11 Comparison Of Drag Curves 
For Different Empirical 

6 

Constants ( Re = 0.28 x 10 ) 
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Figure 4,12 



Comparison Of Transition Length 

For Different Empirical Coefficients G 

6 

At Upper Surface ( Re = 0.28 x 10 ) 




83 



GAMMATR GAMMATR GAMMATR 

0.0 1.0 0.0 1.0 0.0 






0.0 0.2 •• 0.4 0.6 0.8 1.0 



x/c 



LEGEND 

-4 

-2 

0 

"2 

4 

6 



LEGEND 

-4 

-2 

0 

”2 

4 

_6 

8 

To ■■■■ 



LEGEND 

-4 



0 

"Z' 

4 

_ 6 _ 

0 

To' 



(c) 




12 



Figure 4.13 



Comparison Of Transition Length 
For Different Empirical Coefficients 
At Lower Surface ( Re = 0.28 x 10^ 







V. CONCLUSIONS 



Cebeci's interactive computer program was applied to 
the Wor tmann-Al thaus FX 63-137 airfoil to show the 
capability of strong v iscous / inv isc id interaction methods to 
predict airfoil flows at low Reynolds numbers. 

From the comparisons with the experimental data, it 
was confirmed that the results are closer to the 
experimental data as the Reynolds number increases. 

Also, much better results were obtained by decreasing 
the empirical constant 

T-tr 

Therefore, it was concluded that the boundary layer 
transition model has an important influence on the 
predictive capability of viscous/ inviscid interaction 
methods . 
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